Cyril NEFZAOUI & Bertrand PATUREL
Présentation des codes et des résultats
Visualisation des données
Imputation par continuité
Modèle linéaire
GAM
GBM
XGBoost
Imputation VS prévision pure
Permutation
Complétion partielle
On récupère les deux jeux de données train (jeu de données d'entrainement) et test (jeu de données de submission).
setwd("C:/Users/patur/OneDrive/Bureau/ENSTA/ENSTA Cours/Projet R/TP")
rm(list = objects())#comme clear all, met a jour des var
train=read.csv("Data/train.csv",header=TRUE) #on charge le tableau d'entrainement train
test=read.csv("Data/test.csv",header=TRUE) #on charge le tableu de résultat test
On regarde tout d'abord a quoi ressemble notre jeu de données Appliances en fonction du temps.
library(ggplot2)
ggplot(train,aes(Posan,Appliances))+geom_point()+xlab("Temps")+ylab("Consommation")
Puis l'on s'interresse au jeu de donnée test. On souhaite comprendre quelle est la plage des relevés de cette table et la comparer aux dates du jeu de données train.
library(xts)
Date_train=as.POSIXct(strptime(train$date,tz="GMT","%Y-%m-%d %H:%M:%S")) #conversion de la date de train
train.xts=xts(train$Appliances,order.by = Date_train) #conversion de la liste train en série temporelle
end=end(train$Posan)[1]
set=which(test$Posan==train$Posan[end])
l=end(set)[1]
first=set[l]
last=end(test$date)[1]
Date_test=as.POSIXct(strptime(test$date,tz="GMT","%Y-%m-%d %H:%M:%S")) #conversion de la date à prévoir de test
num=last-first
c=rep(200,num)
test.xts=xts(c,order.by = Date_test[(first+1):last]) #conversion de la liste à prévoir du test en série temporelle
res.xts=cbind(train.xts,test.xts)
plot(res.xts) #en noir ce que l'on connait en rouge ce que l'on doit prévoir
Conclusion : On suppose donc qu'il y avait un jeu de données complet initialement qui à été répartie en deux jeux de données, l'un train et l'autre test
Avant tout traitement de train et test, nous avons rajouté ces data frame d'une colonne , dont les coefficients valent 1 ou 0 selon que l'élèment associé appartient au train ou au test. On a ensuite fusionné train et test dans une large que l'on appelle total. On effectue un tri de total par ordre chronologique croissant, ce qui nous permet de visualiser la répartion de certaines valeurs du test dans le train.
#### On reconstitue le jeu de données initial
library(mltools)
library(data.table)
#On créer une variable qui vaut 0 si le relevé est issue de train, 1 si c'est dans le test
yes_is_in_test=rep(1, times = length(test$date)) #On créer la varaible de 1
test$is_in_test=yes_is_in_test #On attache cette variable au test
not_in_test=rep(0, times = length(train$date)) #On créer la varaible de 0
train$is_in_test=not_in_test #On attache cette variable au train
#On reconstitue le tableau initial qui a donné forme aux deux jeu de données
train_dt=setDT(train, keep.rownames=FALSE, key=NULL, check.names=FALSE)
test_dt=setDT(test, keep.rownames=FALSE, key=NULL, check.names=FALSE)
total= merge(train_dt,test_dt,all=TRUE) #On recolle nos deux tableaux ensemble
#length(which(is.na(total$Appliances))) les valeurs manquantes des appliances du test à compléter figure comme des NA, il y en bien 5771
#length(total$lights) vaut 19735
#Puis on réeordonne la série totale en fonction des dates
total$date <- as.Date(as.character(total$date),tz="GMT",format="%Y-%m-%d %H:%M:%S")
total=total[order(total$date,total$NSM, decreasing=FALSE),]
total=total[1:(19735-(5771-4654)),] #On prend seulement les dates à imputer
On cherche en priorité à savoir ou sont les données manquantes dans les jeu de données.
#summary(train) Nous donne les colonnes avec les NA, on remarque donc qu'il y a seulement les colonnes RH_6 et Visibility
#summary(test) Nous donne les colonnes avec les NA, on remarque donc qu'il y a seulement les colonnes RH_6 et Visibility
On remarque que les colonnes RH 6 et Visibility qui correspondent à la visibilité et à l'humidité sur le facade extérieure de la maison ont de valeurs manquantes. Affichons ces deux variables pour comprendre pourquoi il y a des données manquantes :
ggplot(data=train) + aes(x=as.numeric(date),y=RH_6)+ geom_point() + xlab('relevés')+ylab('Visibilité')+ggtitle('Mesure de la visibilité en fonction des relevés')
On comprend donc que les valeurs manquantes sont dus à la saturation des capteurs de l'humidité et de la visibilité lorsqu'il y a des conditions extrêmes.
index=which(is.na(train$RH_6))
#mean(train$RH_6[index-1],na.rm=T) vaut 12
#var(train$RH_6[index-1], na.rm=T) variance faible
summary(train$RH_6[index-1])
index=which(is.na(train$Visibility))
#mean(train$Visibility[index-1],na.rm=T) proche de 10
#summary(train$Visibility[index-1])
print(train$Visibility[index+1],na.rm=T,na.print="Saturation")
#train$Visibility[index+1]
On choisit alors deux méthodes pour compléter les valeurs saturées.
sequence1=which(is.na(train$RH_6))#Nettoyage des données de RH
for (i in sequence1)
{train$RH_6[i]=10} #Pour les NA on prend comme valeur 10
sequence1=which(is.na(test$RH_6))#Nettoyage des données de RH
for (i in sequence1)
{test$RH_6[i]=10} #Pour les NA on prend comme valeur 10
sequence2=which(is.na(train$Visibility))
for (i in sequence2)
{
k=0
j=0
while (is.na(train$Visibility[i+k])==is.na(train$Visibility[i]))
{k=k+1}
while (is.na(train$Visibility[i-j])==is.na(train$Visibility[i]))
{j=j+1}
train$Visibility[i]=0.5*(train$Visibility[i-j]+train$Visibility[i+k])
}#On prend comme valeur pour les NA la moyenne des bords
sequence2=which(is.na(test$Visibility))
for (i in sequence2)
{
k=0
j=0
while (is.na(test$Visibility[i+k])==is.na(test$Visibility[i]))
{k=k+1}
while (is.na(test$Visibility[i-j])==is.na(test$Visibility[i]))
{j=j+1}
test$Visibility[i]=(k/(j+k))*(test$Visibility[i-j])+(j/(j+k))*(test$Visibility[i+k])
}#On prend comme valeur pour les NA la moyenne des bords
Nous venons de mettre en évidence qu'il y avait deux parties, une partie comprenant l'imputation statistique et l'autre de la prédiction. Dans cette partie nous nous concentrons sur une méthode concise.
Notons la série de toutes les Appliances, celle du test à compléter et celle du train. Pour chaque valeur Appliances manquante dans le test, on regarde les deux valeurs de train et encadrant et on effectue une moyenne pondérée de ces deux valeurs. Puis on améliore cette aproximation en attribuant un poids différent à ces deux Appliances encadrantes. Ainsi plus la valeur manquante à compléter du test est proche temporellement de son bord gauche (resp droit) plus le poids attribuer au bord gauche (resp droit) est important.
On utilise donc la relation suivante : =
library(ggplot2)
res=test$rv1
p=1 #Initialisation
sequence3=which(is.na(total$Appliances))
for (i in sequence3)
{
k=1
j=1
while (is.na(total$Appliances[i+k])==is.na(total$Appliances[i]))
{k=k+1}
while (is.na(total$Appliances[i-j])==is.na(total$Appliances[i]))
{j=j+1}
total$Appliances[i]=(k/(k+j))*(total$Appliances[i-j])+(j/(j+k))*(total$Appliances[i+k])
res[p]=(total$Appliances[i])
p=p+1
}
ggplot(data = total)+geom_point(size=0.01)+aes(x = date,y = Appliances)
library(MASS)
level=13000 #On choisit comme donnée d'entrainement les 13000 premières valeurs du train
# On décide de retirer les variables redondantes et non numérique qui n'auraient pas de sens dans une regression
Data=train[1:level,-1] #On enlève la date qui ne peut être lu
Data=Data[,-33]
Data=Data[,-30]
Data=Data[,-30]
Data=Data[,-30]
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)
#Trois manière de trouver le modèle linéaire minimisant l'AIC
step.modelf=stepAIC(fitstart,direction="forward",scope=formula(fitall),trace=FALSE)
step.modelb=stepAIC(fitall, direction = "backward",trace=FALSE,steps=50)
step.model=stepAIC(fitall, direction = "both",trace=FALSE,steps=50)
#step.modelf$coefficients
#step.modelb$coefficients
#step.model$coefficients
# On décide de retirer les variables redondantes et non numérique qui n'auraient pas de sens dans une regression
Databis=train[level:13964,-1] #On enlève la date qui ne peut être lu
Databis=Databis[,-33]
Databis=Databis[,-30]
Databis=Databis[,-30]
Databis=Databis[,-30]
#Databis=Databis[,-1] #On enlève les appliances
step.model.forecastf=predict(step.modelf, newdata=Databis)
step.model.forecastb=predict(step.modelb, newdata=Databis)
step.model.forecast=predict(step.model, newdata=Databis)
source("function/rmse.R")
#On peut comparer les trois rmse en les méthodes de recherche
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf) : 80.8521408523872
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastb) : 80.8443308079354
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecast) : 80.8443308079354
plot(train$Appliances[level:13964], type='l',ylab='original',xlab='temps')
lines(step.model$fitted.values, col='red')
On cherche désormais à ajuster le nombre de degrés de libertés du modèle de la régression linéaire. On se propose tout d'abord un petit code pour voir le nombre de degrés de liberté minimisant le rmse sur un échantillon proche des valeurs à prédire.
rmsef=rep(0,40)#creation d'un vecteur initial de
for (i in 0:39)
{
step.modelf=stepAIC(fitstart, direction = "forward",scope=formula(fitall),trace=FALSE,steps=i)
step.model.forecastf=predict(step.modelf, newdata=Databis)
rmsef[i+1]=rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf)
}
c=seq(0,40)
plot(rmsef,type='l')
On cherche donc à calculer la trace de la matrice . On obtient alors 12 degrés de libertés, l'on choisit donc de prendre 12 varaiables explicatives.
source("function/rmse.R")
X=as.matrix(Data[-1])
M=X%*%(ginv(t(X)%*%X))%*%t(X) #On cherche à connaitre le nombre de degré de liberté
deg=sum(diag(M))#le nombre de degrés de libertés est 12
step.modelf=stepAIC(fitstart, direction = "forward",scope=formula(fitall),trace=F,steps=deg)
step.model.forecastf=predict(step.modelf, newdata=Databis)
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf) le rmse est de 78.9374621296323
plot(train$Appliances[level:13964],type='l')
lines(step.model.forecast,col='red')
On remarque que les données ont du mal à absorber les pics. La regression linéaire à du mal à prédire les valeurs extrêmes.
Remarque : On voit clairement que notre prédiction peut être améliorée
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)
step.model <- stepAIC(fitstart,scope=formula(fitall), direction = "forward",data=Data,steps=12,trace=F)
step.model.forecast <- predict(step.model, newdata=test[(first+1):5771,])
res[(first+1):5771]=step.model.forecast
ggplot()+geom_line(aes(as.numeric(test$date),res))
Finalement on a un score Kaggle - RMSE : 67.9
Une fois notre meilleur modèle linéaire choisit on décide de l'améliorer en classifiant les données. Nous montrons plus tard l'importance de certaines données, notamment la température de la cuisine et son humidité (qui indiquent clairement que la famille est en train de préparer le repas par exemple) où la température et l'humidité de la salle de bain. Si ces varaiables ont des valeurs fortes nous sommes sur des pics pour appliances et ainsi nous boostons les pics de notre régression linéaire à l'aide d'une méthode proche de celle des K plus proche voisins.
Notre code comporte 4 étapes :
data=total[15000:18000,]
#La cuisine est elle utilisée ?
ind=which(data$T1>25.6)
ind2=which(test$T1[4655:5771]>25.55)
pic4=mean(data$Appliances[ind])
res[4655:5771][ind2][5]=pic4
#La douche ou baignoire est elle utilisée ?
ind=which(data$RH_5>70 & data$T5>25.1)
ind2=which(test$RH_5[4655:5771]>70 & test$T5[4655:5771]>25)
pic5=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic5
#Est ce qu'un habitant fait du repassage ?
ind=which(data$RH_7>40 & data$T7>25.1 & data$NSM>=63000 & data$NSM<=63600)
ind2=which(test$RH_7[4655:5771]>40 & test$T7[4655:5771]>24.5)
pic7=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic7
### Si il fait beau le weekend ils vont dehors
ind=which(data$lights==0 & data$WeekStatus=='Weekend' & data$T_out>23)
ind2=which(test[4655:5771,]$lights==0 & test[4655:5771,]$WeekStatus=='Weekend' & test[4655:5771,]$T_out>23)
picw=mean(data$Appliances[ind])
res[4655:5771][ind2]=picw
#Les gens dorment ils ?
ind=which(data$NSM>84600)
ind2=which(0<test[4655:5771,]$NSM & test[4655:5771,]$NSM<3600)
pic1=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic1
ind=which(0<data$NSM & data$NSM<3600)
ind2=which(0<test[4655:5771,]$NSM & test[4655:5771,]$NSM<3600)
pic1=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic1
ind=which(3600<=data$NSM & data$NSM<14000)
ind2=which(3600<=test$NSM[4655:5771] & test$NSM[4655:5771]<14000)
pic2=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic2
ind=which(14000<=data$NSM & data$NSM<24600)
ind2=which(14000<=test$NSM[4655:5771] & test$NSM[4655:5771]<24600)
pic3=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic3
ggplot()+geom_point(aes(x=as.numeric(test$date[4655:5771]),y=step.model.forecast),size=0.1)+geom_line(aes(as.numeric(x=test$date[4655:5771]),y=res[4655:5771]))+xlab('date')+ylab('Prédiction')
Finalement on a un score Kaggle - RMSE : 67.8
library(corrplot)
library(ellipse)
library(qgraph)
corr=data.frame(Appliances=train$Appliances,NSM=train$NSM,lights=train$lights,T1=train$T1,RH_1=train$RH_1,RH_5=train$RH_5,T5=train$T5)
corr=as.matrix(corr)
qgraph(cor(corr))
On décide alors de choisir en partie NSM et lights comme variables explicatives (bien que deux variables corrélés n'implique que l'une explique l'autre) pour commencer. Nous avons ensuite complété notre choix de variables avec un StepGam. Voici en noir la courbe reélle d'une partie des valeurs de test ayant servi d'entrainement et en rouge la prédiction de ses valeurs par le gam.
### Step.Gam ###
library(gam)
Gam.object=gam(step.modelf,data=train)
step.object <-step.Gam(Gam.object,trace=F,direction='forward',scope=list("lights"=~1+lights+s(lights),"RH_1"=~1+RH_1+s(RH_1),"RH_2"=~1+RH_2+s(RH_2),"RH_3"=~1+s(RH_3),"RH_4"=~1+s(RH_4),"RH_5"=~1+s(RH_5),"RH_6"=~1+s(RH_6),"RH_7"=~1+s(RH_7),"RH_8"=~1+s(RH_8),"RH_9"=~1+s(RH_9),"RH_out"=~1+RH_out+s(RH_out),"Visibility"=~1+Visibility+s(Visibility),"Posan"=~1+Posan+s(Posan),"NSM"=~1+NSM+s(NSM),"T1"=~1+T1+s(T1),"T2"=~1+T2+s(T2),"T6"=~1+T6+s(T6),"T7"=~1+T7+s(T7),"T_out"=~1+T_out+s(T_out),"BE_load_actual_entsoe_transparency"=~1+BE_load_actual_entsoe_transparency+s(BE_load_actual_entsoe_transparency),"BE_load_forecast_entsoe_transparency"=~1+BE_load_forecast_entsoe_transparency+s(BE_load_forecast_entsoe_transparency)))
source("function/rmse.R")
predict_gam=predict(step.object,train[10000:13964,])
# rmse(train$Appliances[10000:13964],predict_gam) de 84
plot(train$Appliances[10000:13964],type='l')
lines(predict_gam,col='red')
Finalement on a un score Kaggle - RMSE : 70
Nous avons améliorer nos deux modèles de prédiction GAM et linéaire en faisant une combinaison de ces mmodèles. L'idée est de laisser le modèle GAM s'exprimer lorsqu'il y a des pics importants et lorsque les habitants dorment (pour combler les défaults).
On effectue donc tout d'abord une prédiction avec un modèle gam expliquant les Appliances uniquement par NSM. On souhaite avoir un modèle prédictif valant zéro lorsque les habitants dorment et valant 1 lorsqu'il y a un pic de consommation. On pose donc .
Nos Appliances prédites seront alors telque :
où correspond à notre meilleur prédiction linéaire des Appliances manquantes du test.
Justification de ce modèle : Lorsque les habitants dorment, est proche de donc équivaut à ce qui est en moyenne la valeur des Appliances de minuit à 7h00 du matin. Quand il y a un pic est proche de 1 et donc correspon à la prédiction linéaire .
#On nettoye les données pour la régression linéaire
Data=train[,-1]
Data=Data[,-33]
Data=Data[,-30]
Data=Data[,-30]
Data=Data[,-30]
#On effectue la regression linéaire
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)
step.model <- stepAIC(fitstart,scope=formula(fitall), direction = "forward",data=Data,steps=12,trace=F)
step.model.forecast <- predict(step.model, newdata=test[4655:5771,])
res[4655:5771]=step.model.forecast
model=gam(Appliances~s(NSM),data=train)
p=predict(model,newdata=test[4655:5771,])
p=(p-min(p))/(max(p)-min(p))
#On effectue le modèle mixte
res2=rep(1,(5771-(4655)))
res2=51.2*(1-p)+res[(first+1):5771]*p
mix=data.frame(Appliances=res2)
ggplot()+geom_point(aes(as.numeric(test$date[4655:5771]),step.model.forecast),size=0.1)+geom_line(aes(as.numeric(test$date[4655:5771]),mix$Appliances))+xlab('date')+ylab('Prédiction')
Finalement on a un score Kaggle - RMSE : 67.9
On peut visualiser un indicateurs de la performance : l'erreur Out Of Bag qui doit normalement converger vers 0 quand on augmente le nombre d'itération dans l'algorithme GBM.
Definition : L'Erreur Of Bag est ou est la prédiction en fonction des variables explicatives et la variable cible (ici Appliances). L'Erreur OOB représente l'écart moyen entre la prédiction et la variable cible, en régression on choisit naturellement pour formule avec .
#On réduit le jeu de données, on garde celle qui nous interresse et on enleve les redondances :
library(dplyr)
train_gbm=(train %>% select(-c(date,Day_of_week,DayType,lights,Visibility,InstantF,Posan,Month,Heure,BE_load_actual_entsoe_transparency,
BE_load_forecast_entsoe_transparency,BE_wind_onshore_capacity,BE_wind_onshore_capacity,BE_wind_onshore_generation_actual,BE_wind_onshore_profile,WeekStatus)))
library(gbm)
object_gbm=gbm(formula =my_formula ,
data =train_gbm, distribution = "gaussian",
, var.monotone = NULL, n.trees = 100,
interaction.depth = 1, n.minobsinnode = 10, shrinkage = 0.05,
bag.fraction = 0.5, train.fraction = 0.5, cv.folds = 20,
keep.data = TRUE, verbose = TRUE, n.cores = 4)
On remarque bien que certaines variables ont une influence relative faible. En appliquant le principe de parcimonie, nous avons donc décidé de retirer ces variables pour ne garder que les plus significatives.
Nous remarquons une similitude : ce sont les mêmes variables qui sont retenues dans le cadre de GBM et dans le cadre du modèle linéaire obtenu par step.AIC. On peut donc entrainer un modèle linéaire obtenu par step.AIC à l'aide d'un GBM.
Nous avons utilisé XGBoost en tant que GBM amélioré.En effet,XGBoost tout comme GBM utilisent le principe de gradient Boosting, mais il nous a semblé que XGBoost proposait une implémentation plus performante que GBM, grâce à notamment:
parallélisation du calcul
tous les paramètres clés de l'apprentissages sont modifiables
utilisation de matrices pleines pour accélerer le temps de calcul.
La modélisation a été concue de telle sorte à limiter l'overfitting
En fait, nous avons utilisé XGBoost comme une interface, qui regroupe en fait plusieurs modèles. Premièrement, on définit une structure de contrôle qui permet de choisir la méthode de validation : nous avons successivement choisi la validation croisée (cv) et l'Erreur Out of Bag (OOB).
Nous nous sommes alors demandé comment les arbres de régression développés dans la phase d'apprentissage allaient mesurer l'homogénéité dans les données de Appliances, variable quantitative.
##importation des différentes bibliothèques
library(tidyverse)
library(lubridate)
library(xts)
library(dygraphs)
library(ranger)
library(dplyr)
library(mltools)
library(data.table)
library(xgboost)
library(caret)
p2=read.csv(file="submission_73.2.csv", sep=",", dec='.')$Appliances[4655:5771]
p3=read.csv(file="submission_73.3.csv", sep=",", dec='.')$Appliances[4655:5771]
p4=read.csv(file="submission_73.9.csv", sep=",", dec='.')$Appliances[4655:5771]
p_moy=(p2+p3+p4)/3
p_moy_original=p_moy
p_moy=list(p_moy)
###Méthode de complétion partielle pour conserver les valeurs d'Appliances déjà prédites et vraisemblables
p_moy_df=setDT(p_moy, keep.rownames=FALSE, key=NULL, check.names=FALSE)
p_moy$to_use=rep(0, times = 1117)
p_moy$to_use[which(p_moy$V1<80)]=1
p_moy$V1[which(p_moy$to_use==0)]=NA
par(mfrow=c(2,2))
plot(ts(p_moy_original))
lines(ts(p_moy$V1),
col="red")
###On cree une copie de total et on l'appelle heavy pour insister sur le fait qu'il contient toutes les variables explicatives.
total_heavy=total
total_heavy_matrix=data.matrix(total_heavy, rownames.force = NA)
total_heavy_numeric= (lapply(total_heavy, function(x) as.numeric(x)))
total_heavy_numeric_df=setDT(total_heavy_numeric, keep.rownames=FALSE, key=NULL, check.names=FALSE)
library(doParallel)
library(missForest)
registerDoParallel(cores=4)
###on creee une variable total_light qui contient uniquement les variables pertinentes.
total_light=(total_heavy_numeric_df %>% select(-c(rv1,rv2,lights,Tdewpoint,RH_6,Visibility,Instant,DayType,BE_load_actual_entsoe_transparency,
BE_load_forecast_entsoe_transparency,BE_wind_onshore_capacity,BE_wind_onshore_capacity,BE_wind_onshore_generation_actual,BE_wind_onshore_profile,Month,Windspeed,Press_mm_hg,InstantF,T2,T3,RH_2,RH_8,Day_of_week,Posan)))
permutation=sample(1:19735)
total_light_permuted=total_light[permutation]
total_with_imputation=missForest(total_light, maxiter = 10 ,ntree = 10, variablewise = TRUE,
decreasing = FALSE, verbose = TRUE,
mtry =sqrt(ncol(total_light)),
replace = TRUE,
classwt = NULL, parallelize = 'forests')
##On recupere les donnees avec les Appliances imputees
data_complete=total_with_imputation$ximp
#######Approche Xgboost
library(xgboost)
#####Définition de la structure de controle
xgb_trcontrol = trainControl(method = "cv", number = 1000,allowParallel = TRUE, verboseIter = TRUE, returnData = FALSE)
#####definitiion du grid
xgbGrid <- expand.grid(nrounds =1, max_depth =50,colsample_bytree = 0.9, eta = 0.1,gamma=0,min_child_weight = 1,subsample = 0.5)
total_light$is_in_test[which(is.na(total_light$Appliances))]=2
to_train=total_light[which(total_light$is_in_test<2)]
y_train=total_light$Appliances[which(total_light$is_in_test<2)]
total_light_wt_App=(to_train %>% select(-Appliances))
X_train=(xgb.DMatrix(as.matrix(total_light_wt_App ),label=y_train))
xgb_model=xgb.train(params=list(booster="gbtree"),data=X_train, trControl = xgb_trcontrol, tuneGrid = xgbGrid, method = "xgbTree",nrounds=100000)
rm(X_test)
X_test=total_light[which(total_light$is_in_test==2)]
X_test=(X_test %>% select (-Appliances))
X_test=xgb.DMatrix(as.matrix(X_test))
predictions=predict(xgb_model,newdata=X_test)
plot(ts(predictions))
reconstitution_test=test_total$Appliances
length(which(is.na(reconstitution_test)))
tokaggle=c(appliances_imputed,predictions)
plot(tokaggle)
total_light$Appliances[which(is.na(total_light$Appliances))]=predictions
length(which(is.na(test_total$Appliances)))
test_total$Appliances[which(is.na(test_total$Appliances))]=predictions
which(is.na(test_total$Appliances))
plot(test_total$Appliances[4655:5771])
which(is.na(total_light$Appliances))
tokaggle=test_total$Appliances[4655:5771]
plot(to_g)
On remarque que le test a été echantilloné de manière linéaire dans le train. On dispose donc d'un échantillonage de bonne qualité, car les modèles peuvent être entrainés sur le train, ce qui aurait été plus difficile si les valeurs manquantes de Appliances avaient été réparties sans homogénéité.
Voici le résultat de l'imputation des données manquantes à l'aide d'une méthode de prévision à l'aide de Xgboost en utilisant une méthode de validation croisée.
Par ailleurs, voici le résultat de l'imputation des données manquantes à l'aide d'une méthode de complétion par forêts aléatoires. Plus précisement, nous avons utilisé le . Dans cette méthode d'imputation, l'erreur Out of Bag converge vers 0, ce qui indique le validité de notre imputation.
Finalement on a un score Kaggle - RMSE : 69.1
Remarque : Nous choisissons de travailler sur les jeu de données entier : on considère que les données imputées font partie des données initiales, donc on peut les réutiliser pour entrainer nos modèles. Nous allons discuter de ce choix dans la suite de notre étude.
Précedemment, nous avons établi que les données imputées étaient très proches des vraies données au sens du RMSE.Dès lors, on fait l'hypothèse que ces données sont plausibles.On peut donc considérer ces données manquantes avec deux points de vue:
1) Ces données manquantes peuvent être complétées par imputation (ACP par missMDA, forêts aléatoires par missForest)
2) Ces données peuvent être prédites à l'aide d'un modèle qui a été entrainé sans prendre en compte ces données.
Sur cette figure, on a représenté 3 séries mises en jeu dans le test
en vert, les valeurs imputées par missForest,
en violet, la partie présente dans le train,
en bleu,la partie de prédiction pure.
Conclusion : Pour améliorer notre score, nous décidons d'améliorer notre prédiction pure dans notre algorithme.
Nous représentons ici plusieurs prédictions, en utilisant plusieurs jeux de paramètres dans l'algorithme XGBoost
De plus, puisque tous les prédicteurs sont numériques,on se ramène en fait à une régression constante par morceaux sur des pavés de ,ces pavés étant déterminés par dichotomie succesives lors de la création de l'arbre de régression
Nous avons constaté que, une fois le problème d'imputation résolu, le problème de prévision des valeurs dans le futur présente une difficulté supplémentaire: on doit prédire un un ensemble continu dans le temps de Appliances ,alors que dans le problème d'imputation, on pouvait compléter les valeurs manquantes grâce aux valeurs des voisins des Appliances manquantes, par exemple par la méthode des plus proches voisins.
Nous avons donc appliqué une permutation sur l'ensemble d'apprentissage afin d'homogénéiser la répartition des valeurs à prédire de facon uniforme.Ainsi, on transforme un problème de prédiction en un problème d'imputation.
En vert est représentée une prédiction réalisée à l'aide d'une régression linéaire classique.
En bleu est représentée une prédiction réalisée à l'aide de notre méthode de permutation-imputation.
Désormais, on a trois possibilités de prédire les valeurs manquantes:
1) Prévision à l'aide d'un modèle (régression, GBM, Xgboost, Random Forest)
2) Imputation, sans permutation.
3) Imputation, avec permutation.
Interprétation : Sur ce graphe, on superpose les résultats issus d'une même méthode d'imputation, avec et sans permutation. Avec permutation, les valeurs à prédire sont dispersées au sein du jeu de données et le lien temporel entre chaque relevés successif à predire est moindre. Ainsi la variance est plus élevée.
L'ACF de la série des Appliances d'entrainement montre bien une corrélation de la série à son passé.
On dispose désormais de prévisions de bonne qualité, c'est-à-dire avec un faible RMSE. On choisit nos 3 meilleures prédictions au sens du RMSE, puis on crée la série qui correspond à la moyenne de ces prévisions. Le signal résultant est intéressant car on retrouve de faible appliances la nuit. Ainsi les valeurs d'Appliances très faibles (seuil Appliances <80) paraissent vraisemblables. Nous choisissons donc de supposer ces valeurs comme vraies. On réduit ainsi le cardinal de l'ensemble des données à prédire.En fixant ce seuil , on ne prédit par exemple que 893 valeurs au lieu de 1117. En rouge sont représentés les valeurs de Appliances qu'on introduit dans le jeu d'entrainement, et qu'on a donc plus besoin de prédire.
On peut comparer les avantages et les inconvenients de nos diverses méthodes. Dans un premier temps, nous avons mis en place une simple regression linéaire.Cette méthode a le mérite d'être facilement interprétable, et elle constitue une méthode qui conserve une continuité.En effet, de petites variations sur les données conduiront à de petites variations sur les résultats, grâce à la continuité du produit matriciel. Cependant, les régressions linéaires sont par nature inadaptées pour représenter des variations brusques (à savoir les pics).
Au contraire, les méthodes d'arbre permettent de modéliser de facon plus fidèle les discontinuités.Cependant, de faibles variation sur les données peuvent conduire à de grandes variations sur les résultats.Cependant, ces défauts ont été écartés grâce à l'utilisation du bootstrap.On a aussi établi que notre cas, l'utilisation d'un arbre de régression s'identifiait à une régréssion linéaire sur des fonctions constantes par morceaux, sur des pavés de .
Remarque générale : Lors de l'utilisation de nos modèles, nous avons été confrontés au compromis biais-variance. Des modèles simples, avec peu de variables explicatives ont eu généralement un fort biais mais une variance peu élevée, et les modèles complets ont parfois réussi à obtenir un biais inférieur mais une variance élevée.
Nous sommes deuxième de la compétition.